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Abstract 

. We analyze the collective behavior of a lattice model of pulse-coupled oscillators. 

By studying the intrinsic dynamics of each member of the population and their 
mutual interactions we observe the emergence of either spatio-temporal structures 
or synchronized regimes. We perform a linear stability analysis of these structures. 
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1 Introduction 
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^ ■ The collective behavior of large assemblies of pulse-coupled oscillators has 

been investigated quite often in the last years. Many physical and biologi- 
cal systems can be described in terms of populations of units that evolve in 
time according to a certain intrinsic dynamics and interact when they reach a 
threshold value [1] . Although it was known long time ago that the members of 
these systems tend to have a synchronous temporal activity, a rigorous treat- 
ment of the problem has been considered only in the last decade [2-4]. Up to 
now, the most important efforts have been focused on systems with long-range 
interactions because in this case analytical results can be derived by applying 
a mean-field formalism. Relevant is the work by Mirollo and Strogatz (MS) 
[4] who discovered under which conditions mutual synchronization emerges as 
the stationary configuration of the population. Later on, the study has been 
generalized to other situations [5-10]. 

When the oscillators form a finite dimensional lattice where only short-range 
interactions are allowed, the spectrum of behaviors is broader. For example, 
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under certain conditions lattice models of pulse coupled-oscillators display self- 
organized criticality (SOC) [11-16]. In such case the system self-organizes, due 
to its own dynamics, into a critical state with no characteristic time or length 
scales and provide some hints about the wide ocurrence of 1/ f noise in nature 
[17]. In other cases, phase locking, synchronization, or more complex spatio- 
temporal structures are developed. 

On the other hand, some efforts have been devoted recently to analyze the 
stability of spatio-temporal periodic structures in coupled map lattices [18]. 
These systems are easier to deal with, both analytically and numerically, than 
lattices of coupled nonlinear oscillators although they represent a less realistic 
view of coupled systems that appear quite often in nature. 

It is precisely the purpose of the present work to study the stability of certain 
structures, that have been obtained recently in computer simulations in a 
lattice of pulse-coupled oscillators. Keeping this goal in mind the paper is 
organized as follows. In the next section we introduce the model and explain 
several equivalent ways to describe the time evolution of the system. These 
equivalent descriptions are related through well defined transformations, which 
allow to study the model in terms of the most suitable dynamic variables. In 
Sect. 3, we calculate analytically the fixed points of some of the structures 
observed in simulations in one- and two-dimensional lattices, and in Sect. 4 
we analyze the stability of these fixed points. Finally, Sect. 5 is devoted to the 
conclusions of this work and its future perspectives. 



2 The model 

Let us consider a system of coupled oscillators described each one by a state 
variable Ei, which we can identify with a voltagelike magnitude when dealing 
with biological oscillators. We will assume that all the oscillators are identical 
and evolve in time according to the following dynamics 



plus the reset condition for E^ > 1. This means that when the i-th oscillator 
fires, it is reset (Ei > 1 — > 0). The first terms of the r.h.s. is the driving rate 
that throughout the paper will be considered a positive quantity, f(E) > 0. 
The second term accounts for the coupling, given in terms of a nonlinear 
coupling function e(E) that we will assume either excitatory for all E, e(E) > 
0, or inhibitory, e(E) < 0, WE (except at the reset point E = 0, where it could 




dEi 



(1) 



be e(0) = 0). 
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When some of its neighbors (labeled by index j) fire at time tj, the i-th oscil- 
lator suffers an instantaneous perturbation of its state given by the coupling 
function e(Ei). This implies the existence of two time scales, a slow one for 
the driving and a fast one for the coupling. In this description both time scales 
appear in the same equation. This represents a quite general evolution equa- 
tion, because not only the driving rate f(E) but also the coupling e(E) are 
functions of the state E. 

We could also consider another two equivalent descriptions of the same model, 
both related to (1) by simple changes of variables. The first equivalence is 
obtained by applying the following nonlinear transformation 



used before in different contexts by several authors [8,19], where Eq is defined 
to ensure that y(E — 1) = 1. By substituting into (1) we get 



Now, the evolution of the system is described in terms of a new variable y for 
which the coupling is constant and whose driving rate is given by 



A case of particular interest is that of a zero advance at the reset point (e(E = 
0) = 0). In such case, this transformation is well defined if the coupling is 
constant (eq) \/y ^ except for y = where is exactly zero. This condition 
plays the role of a refractory time in the fast time scale, which provokes that 
all the units which have fired have zero phase at the end of the interactive 
process. 

The second equivalent description can be derived by writing (1) as a function 
not of the state Ei but of the phase of each oscillator which evolves linearly 
in time except when it receives a pulse from its neighbors. In this description, 
the driving and the coupling are integrated in a 'phase response curve' (PRC), 
function that gives the phase advance of the oscillator that receives the pulse 
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(3) 




(4) 



[20,21], 
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This equation can be obtained in a straightforward manner by applying the 
following transformation in (1) 



E 



o 



f(E' 



The fact that <j)(E = 1) = 1 is guaranteed if the intrinsic period of the os- 
cillators is equal to one. According to the definitions the function A is given 
by 



All three descriptions are physically equivalent, but depending on the the- 
oretical framework it is convenient to deal with one or another for specific 
purposes. In this paper we will use the description in terms of the phase. Fi- 
nally, it is worth noting that one has to be very careful when handling with 
the two coexisting time scales in the system. Since they do not overlap the 
situation is equivalent to stop the driving when an oscillator fires. Then, the 
resulting PRC A n (0) must be related with the function A(0) by 



4>+&n{4>) tj 



A(0') , 

«7 



where we obtain a different PRC depending on the number of firings n that 
the oscillator receives, which will be bounded between one and the number 
of neighbors. This ensures that when n neighbors fire simultaneously, the i-ih 
oscillator modifies its phase as 

0. + A n (&)- (9) 



otherwise it evolves as dtfii/dt — 1. In general, to go from A(0) to A n (</>) is 
straightforward, but the inverse implies to solve an integral equation that in 
general can only be done numerically. 



3 Fixed points 



First of all, we will show that a lattice model of pulse-coupled nonlinear oscil- 
lators evolving in time following (9) with nearest-neighbor interaction and pe- 
riodic boundary conditions has some periodic spatio-temporal solutions which 
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correspond to fixed points in a return map description. With this goal in 
mind, we will study initially the behavior of a system of two oscillators and 
after this we will generalize the analysis to one- and two-dimensional lattice 
models. Hereafter we will consider always refractory time in the model, i.e., 
A n (0 = 0) = 0. 



3.1 Two oscillators 



To start the discussion we will consider a system of two oscillators. It is simple 
enough to allow a complete description of the dynamical evolution of the model 
and, additionally, to illustrate the main ideas that we want to develop later 
on. We present here a detailed evolution of the phase of each oscillator for half 
cycle. Each row corresponds to the phase of each one. 

1 — ► — ► 1-0-AxO) 

firing driving (10) 

— > + — > 1 



The fixed points of this transformation are the solutions of 0* = 1 — 4>* — 
Ai(0*). If Ai(0) is a continuous function bounded between -1 and 1, there 
exists at least one fixed point. However an important question arises about the 
uniqueness or multiplicity of those fixed points. It is easy to know under which 
sufficient conditions the fixed point is unique by simple geometrical arguments, 
which is A'^0) > —2, V0. It can be proved, with our hypothesis of f(E) > 
and e(E) > or e(E) < 0, WE, that A' n (<j>) > -1, V0, verifying thus the 
uniqueness of the fixed point. To clarify this point, let us consider a particular 
example that allows to calculate the location of the fixed point, the Peskin's 
model. It was proposed to analyze the collective behavior of an assembly of 
cardiac pacemaker cells [22]. In this model the intrinsic time evolution of each 
member of the population is given by 

f(E)= 1 (K-E), (11) 



where 7 gives the slope of the driving rate and K = (1 — e 7 ) . The coupling 
function is for this simple case a constant e(E) = e . Then the function A is 

^ = W^E) - ^ < 12 > 



where for the last equality it is necessary to know the relation E((f>), substi- 
tuting Eq. (11) into Eq. (6). Now it is easy to find the PRC according to (8), 
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obtaining 

A n (0) = -iln(l-^e^), (13) 

It is straightforward to see from the derivative of A n (0) that, for a given 
Eo > 0, the PRC is an increasing function when 7 > 0, whereas for 7 < the 
PRC decreases monotonously. This behavior is the opposite when an inhibitory 
coupling is taken into account. For two oscillators the fixed point is unique 
and satisfies Ai(0*) = 1 — 20*, corresponding to 

^^^-^sinh-^eosinhQ)). (14) 

In general for the two oscillators case we can study the global stability of the 
system [4]. To show this, we look at the result obtained after half a cycle 
assuming identical oscillators. Then for A' x (0) > 0, any perturbation 5 from 
the fixed point evolves according to (10) as 

o = 0* + 5 — >i-0*-5- + <(/)*- 5 

S v ' 

>Ai(4>*)=l-2</>* 

and 

0o = 0* - S — » 0* + 5. 

This behavior shows that the fixed point is a repeller of the dynamics, when the 
PRC is an increasing function, but due to the periodicity of the variable this 
fact provokes synchronization between oscillators, i.e., oscillators are repelled 
from the fixed point and they go to phase one or zero which are the same cyclic 
point. The two oscillators approach each other after each cycle until they are 
close enough so that the firing of one of the oscillators makes the phase of the 
other one to reach the threshold. At this time, due to the refractory time, both 
oscillators fire in unison and this situation persists forever, which corresponds 
to the synchronized state. In the same way, if A' : (0) < we have: 

o = 0* + 5^ >0*-5 
0o = 0* -<f— > <<j)* + 5 

Both equations need another limit 

0* + 5^1-0*-5-A!(0* + 5) <0* + 5 
0* - S -> 1 - 0* + 6 - Ai(V - 8) > 0* - 8 
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where we have use that Ai(0) > 1 — 20 if > 0* and Ai(0) < 1 — 20 if 
< 0*. Now, for a decreasing PRC, the fixed point is an attractor of the 
dynamics. This attractor corresponds to a phase-locked state that maintains 
the difference of phases between both oscillators. Note that we have assumed 
that the initial conditions are chosen in such a way that in the first interaction 
the oscillators do not synchronize. This kind of synchronization is trivial and 
does not depend on the slope of the PRC, only on its strength. 

We have to mention, to clarify this point, that a PRC is always bounded by 
A n (0) < 1 — when it is excitatory and by A n (0) > —0 when inhibitory, in 
order to avoid phases larger than 1 or smaller than 0. However, these limits 
correspond to trivial synchronization of the oscillators, and will not be taken 
into account here. Therefore, when we say increasing or decreasing PRC we 
refer to the domain in which synchronization is not trivial. 



3.2 One- dimensional lattices 

We have seen in the previous paragraphs that depending on the slope of the 
PRC two oscillators tend either to synchronize or to keep a phase difference 
between them. When generalizing this result to a periodic one- dimensional 
model of coupled oscillators with nearest-neighbor interaction, we can imagine 
the same tendency. On the one hand the strength of the coupling makes two 
neighboring oscillators with close phase values to synchronize. On the other 
hand, when they are not close enough, this coupling can make the phases to 
approach each other or to separate, depending on whether the slope of the PRC 
is positive or negative, respectively. We have observed in computer simulations, 
starting from a random distribution of phases, that the population reaches, 
after a long transient and depending on the parameters of the system, well 
defined spatial structures, with concrete values of the phase differences. We 
have checked that these structures are limit cycles for the dynamical evolution 
of the lattice. When building up the return map, by looking to the system at 
the time that a given oscillator reaches the threshold, we obtain fixed points. 

The most simple fixed point is that one in which all the oscillators have exactly 
the same phase. Thus all of them reach the threshold simultaneously. Having 
assumed that A„(0 = 0) = 0, i.e. the existence of a refractory time, all the 
oscillators will be at zero phase after the collective firing and then a new cycle 
starts again. This corresponds to global synchronization of the population. 

There are other fixed points characterized by the fact that the population is 
split into several subpopulations with the same number of elements. Each sub- 
population is composed by oscillators with the same phase, which is different 
for different subpopulations. Moreover, they form periodic spatial structures 
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Fig. 1. General Id spatially periodic structure. Each <pi corresponds to a different 
value of the phase. 



Fig. 2. Id 2-phase fixed point structures: Top) Chessboard-like (2a); bottom) 
Domino- like (lals). 

and have been discussed in the context of coupled map lattices [18]. A general 
case is plotted in Fig. 1, where the subscripts stand for the different subpop- 
ulations. Among these structures we will focus on the simplest ones, those 
having only two different phases, plotted in Fig. 2. We have introduced a no- 
tation to denote the number of synchronized (s) or phase-locked (a) neighbors 
of a given oscillator that will also be used for the 2d models. Without loss 
of generality we can assume that white sites in Fig. 2 have phase equal to 
one whereas black sites have a phase 0. In these cases it is very simple to 
compute the value of 0* that correspond to the fixed points. The proof runs 
parallel to that for two oscillators. For the chessboard, if each row stands for 
a subpopulation 

1 — ► — > l-0-A 2 (0) 

firing driving (15) 

— > + A 2 (0) — > 1 



and the same for the domino writing Ai(0) instead of A 2 (0). Then the fixed 
points correspond to 



chessboard: 0* = 1 - 0* - A 2 (0*) = 1 - 0* - A* 2 
domino: 0* = 1 - 0* - A : (0*) = 1 - 0* - A* 

In the next section we will analyze in detail the linear stability of the chess- 
board structure (2a) in Id and how to generalize this result to other structures. 
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3.3 Two-dimensional lattices 



In two-dimensional lattices, with periodic boundary conditions and nearest- 
neighbor interactions, we have also observed in computer simulations the ex- 
istence of well-defined patterns that are fixed points for the discrete dynamics 
of the population. Some of the two-phases structures are plotted in Fig. 3, to- 
gether with the unique one-phase structure, the synchronized state. The phases 
of the fixed points are very easily computed using the transformations for the 
fixed point explained in the one-dimensional case. In general we are going to 
have for the different possibilities between phase-locked (a) or synchronized 
(s) neighbors in a 2d square lattice, the following fixed point equations: 



4a 


0* = 1 - , 




3als 


0* = 1 - . 


f-A; 


2a2s 


0* = 1 - , 


f-A* 


la3s 


</>* = !-, 


f - A* 



(16) 



In our computer simulations we have observed that these structures have 
different basins of attraction depending on the values of the parameters as 
well as on the size of the lattice. Nevertheless, the chessboard (4a) seems to 
be the most stable in two senses: it has the largest basin of attraction and it 
is the most robust to different sources of dynamical noise. In the next section 
we will briefly discuss the linear stability of these patterns. 



4 Stability of the periodic structures 

Now we will analyze the stability of the fixed points discussed in the previous 
section. We will discuss in detail the Id chessboard, and later on we will explain 
how to generalize this to other structures. 

4-1 Linear stability analysis of the one- dimensional chessboard 

To perform the stability analysis in the one-dimensional model we are going 
to deal with the most simple configuration of two phases (black and white) 
slightly perturbed at random from its fixed point (..., 1, 0*, 1, 0*, ...). Suppose 
without loss of generality we have the subpopulation of white sites distributed 
slightly below phase one, and the subpopulation of black sites with phases 
around 0*, that is 
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Fig. 3. 2d 2-phase fixed point structures and the synchronized state (4s). 
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White sites (W) 



1 first 
I — Si last 

(17) 



Where the top row stands for the first oscillator that will reach the threshold 
and the bottom row for the last one, after a first driving of amount Si. On the 
other hand we have for black sites 



, X + 5 2 first 
Black sites (B) { (18) 

i* + <5 3 last 

where the two rows have the same meaning as for the white sites. Nevertheless, 
this does not mean that these oscillators have the largest and smallest phases, 
in contrast to the white oscillators. Since during the driving Si the black sites 
will receive pulses from their white neighbors at unknown phases, the order 
is not maintained and the black oscillators can overtake each other. Then the 
two oscillators shown in (18) are those that after receiving the two pulses have 
the maximum and minimum phases, respectively. 

As we have stated before, the first step is to wait for a driving Si, at this point 
all the white sites have fired and the new configuration becomes, if in order 
to simplify we call A to Ai: 



W 



B i 



Si 


0* + S 2 + Si + A(0* + S 2 + otiSi) + A(0* + S 2 + a 2 5i + A(0* + S 2 + a^)) 
0* + S 3 + Si + A(0* + S 3 + (3iSi) + A(0* + S 3 + (3 2 Si + A(0* + S 3 + /Mi)) 



where aii, a 2 , j3i, f3 2 are the unknown fractions of the driving Si the oscillators 
have run at the moment they receive the pulses from the white neighboring 
sites. At this point we linearize A around 0*: 

A(0* + 5 2 + aiSi) = A(0*) +(5 2 + a^i) A' (</>*) 

A* A'* 

A(0* + S 2 + a 2 5x + A(0* + S 2 + aiSi)) = 

= A(0* + A*) + (S 2 +a 2 S l + {S 2 +a l S l )A'*) A'(0* + A*) 

S v ' V v ' 

With this linearization, the previous structure can be written for the black 
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oscillators as 



>* + A* + A** + C5 2 + AS 1 / x 

B (19) 
0* + A* + A** + C5 3 + B5 1 



where A, B, and C have the following expressions: 
A = 1 + «iA'*(l + A'") + a 2 A'** 

5 = 1 +/3iA'*(l + A'**) +/3 2 A'** (20) 
C = A'* + 1 + A'** + A'* A'** = (1 + A'*)(l + A'**) 

Now, we drive the population until the first black oscillator arrives to the 
threshold, that will be a driving 1 - (0* + A* + A** + C5 2 + ASi), equal to 
0* — ASi — B5 2 , due to the fixed point condition. The new configuration is 
thus 



B 



0* + (1 - A)8 1 - C5 2 
0* - A5 1 - C5 2 
1 

1 - ((A - B)5 1 + C(5 2 - 5 3 )) 



These steps explain exactly half a cycle of the process. Then we can define a 
new set of perturbations as 



5[ = -C(5 3 -5 2 )-(B-A)5 1 
5' 2 = S\ — C5 2 — A5\ 
$3 = —C5 2 — ASi 



(21) 



Now we can compare with the initial configuration. The white oscillators are 
bounded between 0* + 5' 3 and 0* + 5' 2 , whereas the black ones are between 
1 — 5[ and 1. Then the transformation matrix for the perturbations 5±, 5 2 , 5 3 
is 



A-B C -C 
l-A-C 
-A -C 



(22) 



The eigenvalues of (22) give the information about the linear stability of the 
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fixed points. They are: 
Ai = -C, 



Now we use that A' ((f)) > — 1 always, and then C > 0. The absolute value of 
the eigenvalues will be 



|Ai| =C, 



A — B\ + ^\A-B\ 2 + AC 



(24) 



Due to the fact that < ctj,/^ < 1 (by definition), we have, for a decreasing 
PRC and using Eq. (20), that C < A, B < 1, that is, | A - B\ < 1 - C. This 
yields 

|Ai|=C<l, |A 2 , 3 | < \ (l - C + v/(l - C? + 4C) = 1 (25) 



The fact that the three eigenvalues are less than one, in absolute value, ensures 
the stability of the fixed point, when the PRC is a decreasing function. 

On the other hand, for A' (<f)) > and using again that < ctj, < 1 and 
Eq. (20), we have 1 < A, B < C, which implies \A — B\ < C — 1, and then 
| Ai| = C > 1, | A2,3 1 < C. The fact of having at least one eigenvalue with 
absolute value larger than one ensures that the fixed point is unstable if the 
PRC is an increasing function. 



4-2 Generalization to other structures 



The analysis of the 'domino' structure in Id follows the same ideas developed 
for the chessboard structure if an excitatory coupling is assumed. The only 
difference is that a given unit receives always a pulse from another unit of 
the same color, and a further pulse from a unit of different color. This fact 
introduces new values for A, B, and C in the transformation matrix (22). 
In particular, the value of C is given by C = (1 + A'*). The analysis of the 
transformation matrix of the perturbations shows that all the eigenvalues Aj 
satisfy |A;| < 1 provided A' ((f)) < 0. As this transformation is iterated n times 
(n — > oo) the perturbations will tend to zero, ensuring the linear stability of 
the structure. 

A generalization of the results for the Id lattice to the 2d is also easy to carry 
out and can be applied to any structure shown in Fig. 3. In these cases, the 
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value of C depends on the number of phase-locked neighbors (a) as 

C= (1 + A'*)(l + A'**)... (26) 
„ ' 

a times 

Again, following the same formalism it is easy to prove that the new structures 
satisfy the same stability criteria as those in one dimension. 



5 Conclusions 

In this paper we have studied the collective behavior of a lattice model of pulse- 
coupled nonlinear oscillators. Each single unit has its own intrinsic dynamics 
and interacts instantaneously with its neighbors when it reaches a threshold 
value. At this moment the oscillator is reset. The response of an oscillator at 
the reset point plays a relevant role when considering the dynamical properties 
of the model. The existence of a refractory time is assumed throughout the 
paper. 

The continuous dynamical evolution of each unit can be transformed into a 
discrete one by looking at the system when a fixed oscillator reaches the thresh- 
old. This picture is appropriate to describe the features of spatially periodic 
patterns which are fixed points of the discrete dynamics, and whose linear sta- 
bility has been analyzed. Our main result is that a monotonously decreasing 
(increasing) phase response curve makes these structures to be linearly stable 
(unstable). However, a detailed analysis of their basins of attraction is still 
missing and deserves future attention. 
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